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SUMMARY 


OBJECTIVE 
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RESULTS 
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appropriate  to  the  vlf  communications  band. 
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I  INI  ROUl't  l  ION 


I  his  i  epoi  I  dosei  ihes  .iiul  1 ! ni  s . i  computet  progium  designed  to  handle  pulse  piopaga 
lu'ii  problems  wlu'n  the  piopagalmg  channel  is  the  e.uth-ionosphere  waveguide  and  is 
intended  lor  use  in  tiro  ell  \ 1 1  hands  Inpuls.ite  mode  data  ti  e  .  eigenangles  ami  excitation 
tactorsl  .is  .1  tuneiion  ol  trei|uoncx  .is  determined.  lot  example.  bx  the  waveguide  piogram 
ot  Ret  I  I  lie  nis'vl e  data  .lie  111  lei  polaleJ  using  etihie  splines  lh.it  is,  the  real  ami  imagmai  v 
part s  ol  the  eigenangles  .is  well  .is  the  in.i.cn it tule  .uul  phase  ol  the  exeil.ilion  factor  are 
approximated  In  a  third  -degree  polx  noinial  between  e.ielr  p.111  ol  data  points  Die  polx  110- 
mials  aie  determined  stieh  that  t hex  lit  the  input  data  ami  are  tw  ue  continuouslx  dilteren- 
liahle  in  (he  domain  ol  interest  l  he  pulse  shape  mtegials  tw  hu  ll  are  I  ouuei  Hailstorms!  arc 
calculated  usury:  a  fast  I  ouriet  trail  slot  i  teehnu|tie  \dxantages  and  disadvantages  ol  the 
last  lounei  transloim  in  pulse-shape  studies  hare  been  discussed  In  Sevier.  Bloch  and  |  1\  nn 
1  Ret  .1  lls  niajoi  advantage  is  a  s.i\  mgs  in  computational  time,  w  here.is  a  disadvantage  max 
be  that,  strictlx  .  onlx  periodic  pulse  trams  ma\  he  aualv/ed  I  Inis,  when  a  non-penodie 
pulse  is  considered,  it  must  he  treated  as  a  peitodie  pulse  tram  with  peinul  much  groatoi  than 
the  pulse  w  ultli  in  otdei  to  obtain  adequaie  lesolulion 

\I  the  present  the  projiram  is  designed  to  calculate  onlx  the  vertical  electric  Held  tl  .) 
at  the  ground  loi  a  ground -based  vertical  electric  dipole  source  Otlici  source  ami  receive! 
orientations  ami  altitudes  can  ho  treated  In  extensions  of  the ‘VIIANI  I  "subroutine 
"C  II  VM  l  could  also  he  extended  tv'  allow  tor  lateral  mhomogoncitv  v'l  the  guide  via  WkH 
or  mv'vle  cv'ntersii'ii  metlu'vls  Die  program  was  developed  |'i miarilx  as  a  tv'ol  tor  calculating 
sU'w  -tail  atnu'sphenc  waveforms  ( i.e  .  w  are  shapes  in  the  ell  h.uul  generated  In  lightning  ills 
charged.  Because  anisotropv  v't  the  u’liosphere  is  included  in  calculating  the  input  mode 
parameters,  the  program  i'  particularlv  suited  lo  stmlies  relating  to  geomagnetic  influences 
i'ii  such  wareloinis  It  can  also  he  used  lo  examine  atmospheric  signatures  m  the  \ It  band 
and  tv'  conduct  pcilormaiice  stmlies  on  proposed  01  existing  spre.nl  spectrum  sv  stems  which 
operate  m  the  ell  \  1 1  hands 

Die  mathematical  problem  at  hand  sitnplx  reduces  to  the  calculation  of  a  louriei 
integral  tin  w  Inch  the  integral  is  made  up  v't  a  1 1  ansni  it  ter  spectrum,  recenei  spectrum,  ami 
channel  spectrum,  each  v't  which  w  ill  he  vliscussevl  more  ftillx  m  the  follow  mg  section  Sec¬ 
tions  III  ami  1\  describe  the  progiam,  Sections  \  ami  VI  ci'iitain  I'litput  vlescrtptn'n  ami 
some  sample  results  I  lie  appendix  contains  a  ptogram  listing. 


11.  Ol  VNN1  I  .  RH  1  IVI  R  AND  SOI  R(  I  MODU  S 

In  this  section  equations  used  loi  calculating  a  receive!  output  xxaxcl'oim.  liyp.O. 
related  h’  the  xei deal  electric  held  at  (ho  ground,  will  he  gix on  In  terms  v't  the  great  circle 
range,  p.  the  waveform  is  gixeu  bx  ti  N  -1 1 
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u  here 


idl(iolr(io>h(u>,pt  e'w,du) 


Ul l  (l  )R(ITH(l  ,p)  cj-^^'dF 


(II 


hi l  il  l  iillliol  =  idl(-trl  ). 

source  function  (anipere-m  ll/t 

(2) 

R(  IT  *  rtul  =  r(  JtrlT. 

receiver  function 

(3t 

lid  ,p)  =  in io.pt  =  h( 2 irl- ,p ). 

channel  function 

(4t 

Hie  second  and  third  equalities  in  1  q 
quantity  vi  that 

(It  follow  from  the  requirement  that  Cl(p.t)  be  ; 

i  real 

Idl  (IT  =  |ldl  t-l-TI*.  R(IT  = 

R*(-l  t  and  H(l'.p)  =  1 1  *(— I'.p ) 

(51 

where*  (lie  asterisk  denotes  the  complex  conjugate.  The  receiver,  source  and  channel  lunctions 
are  described  below: 

RECEIVER 


R1  CVRtlT  is  a  subroutine  which  can  be  easily  modified  or  replaced  to  accommodate 
the  individual  users  needs.  I  he  particular  Rl  ('VR(l  l  subroutine  contained  in  the  program 
listing  in  tlie  appendix  can  be  used  with  any  receiver  function  of  the  form 


rue)  = 


(<') 


w  here  the  gain.  CIA.  angular  frequencies  to  |  and  and  integers  P  and  Q  are  read  into  the 
program  via  namelist  This  receiver  function  description  allows  for  a  broad,  but  by  no  means 
exhaustive,  class  of  realistic  receivers.  Observe  that  it  satisfies  the  condition  specified  in 
I  q.  1 5 >. 

I  R  ANSMIT  ft  R 


I  R\M  1'RdT  is  a  subroutine  which  too  can  be  readily  altered  to  meet  the  specific 
needs  of  the  user.  Since  the  principal  motivation  for  the  present  program  was  to  study  the 
shape  of  slow  wave  tails  associated  with  atmospheric  discharges,  the  particular  source  func¬ 
tion  contained  in  the  subroutine  I'RXM  I  Rtl  1  in  the  program  listing  in  the  appendix  is  the 
Williams  (Ref.  3)  mean  source  description  for  a  lightning  discharge,  which  is  given  by 
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1  ho  Aj's,  •>,  s.  fp.  rv  ami  v0  are  ro;ul  into  tho  program  via  namolisl.  llio  units  of  amporos  tor 
1  ho  A  s  atul  m  soo  ti>r  \  0  oouplovl  with  tho  oh  anno  I  slot  moil  in  tho  follow  mji  snhsoot  ion  v  iolils 
a  waveform  m  units  of  volts  in.  vv  liioh.  as  stated.  is  related  to  tho  vortioal  electric  field  at  tho 
mound  I  or  a  Hat  receiver  tho  wavot'orm  wouhl  ho  proportional  tv'  tho  fiolvl,  tho  proportion¬ 
ality  constant  hoinjs  tho  rooonor  gain,  \ 

t'HANNt'L 

1  ho  ohannol  tunotiv'ii  tv'i  tho  oarth-ionosphoro  waveguide  with  a  jiround-hased  vorti¬ 
oal  olootrio  vlipv'lo  sv'iiroo  is 
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o  speed  v't  lijiht  in  vacuum  thin  soo  I 
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,t\\ 


,t\\  \ 

,1(»  /„ 


vortioal  vlij'v'lo 
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w  hero 


l«l 
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llK  II 


1  rosnol  1  M  ground  reflection  ov'olfioiont 

I  rosnol  II  ground  reflection  coefficient 

I I  piano  wave  ionospheric  reflection  coefficient 
1  XI  piano  wave  ionospheric  reflection  coefficient 


ii 


I|K[  I’M  to  II  plane  wave  iouospheric  conversion  coefficient 
[K||  -  II  to  I’M  plane  wave  ionospheric  comer  sion  coefficient 


dW  \ 
•tit 


derivative  of  modal  equation  .1 1  eigen.mgle  0 


"  * 1  ~  II ^  II  11^  ll' ' '  "  1  I*l  1^1 >  _  ||  ||K[  ||R  ||  =  modal  function 

Ml  reflection  ami  am  version  coefficients  arc  referenced  to  tlio  ground.  I  he  channel  func- 
tion.  Ild.pl.  is  defined  sikIi  that  tlio  wave  form  associated  with  the  vortical  electric  field  at 
the  ground 


Up. ti  =  'K^  I  Ml  tl  )Ktl  illd  .pi  e1'*'  1  dl 

Jo 


(ID 


is  in  volts  in  when  hi  It  Id  ulltccl  idlUirl  I  is  in  ampere-in  11/  In  the  lower  elf  hand, 
except  within  a  tew  degrees  ot  broadside,  Fq.  (D)  can  he  used  for  a  ground-based 
horizontal  dipole  radiator,  such  as  the  Wisconsin  Test  Facility,  if  the  excitation  factor 
•or  the  single  mode  (termed  the  N=1  mode)  that  propagates  at  that  frequency  is  replaced 
by 


w  here 


■\  1  cos  p  (  \CS1  |  I 


1  is  given  In  l;q.  1 1 0) 


.  .Of 

is.,  v  / -  *  —  complex  ground  retractive  index 

,u,to  lo 


d:> 


tl.d 


w  here 


(i  -  ground  conductivity  (siemens  nil 

t  =  ground  permittivity  (farads  ml 

t0  =  permittivity  of  free  space  =  S.S5  \  ID-'  -  farads  m 

and  w  here  C  is  the  angle  between  the  direction  ot  the  horizontal  dipole  and  the  direction  of 
propagation 


III  INSCRIPTION  OF  INPUT 

Ml  input  to  the  pulse  shape  program  is  read  in  via  the  card  reader.  \  listing  of  a 
sample  input  showing  the  data  deck  set-up  is  given  on  pages  v>  and  10.  This  sample 
input  applies  to  a  single-mode  case 

There  are  three  parts  to  the  input  I  he  first  part  is  plot  identification.  The  second 
part  is  general  input  read  in  by  means  of  a  namelist  format  The  third  part  is  mode  data. 

I  ach  part  will  be  discussal  in  further  detail  below 


I  lk-  lust  p.iit  ol  Ilk-  input  consists  ol  three  cards  containing  plot  l.thel  information. 

\ll  ear el s  .ire  read  m  using  .1  111  \ 4  torin.it  I  lie  In  si  e.iul  eontains  from  I  to  40  alphanumeric 
eharaeters  containing  vvhatevet  111I01  niation  the  usei  w  ishes  to  he  printed  on  the  transniittet 
speetrum  plot  Ihe  seeoiul  anil  third  card-,  are  iilenlk.il  in  format  ami  eontam  information 
tor  the  receiver  speetrum  plot  anil  the  channel  speetrum  plot,  respect  ivelv 


I  he  seeoiul  part  of  the  input  is  re.ul  in  In  means  of  an  ASCII  I  OK  I  R  AN  namelist 
input  format.  I  he  following  variables  aiul  arrays  mav  he  speeifieil  in  the  namelist  input 


NM 
m  ' : 

M  \  I  1 1  I  I 
\i  vi 4:1  /  1 
M  \  l  (.0  '  I 
M  \  I  'l  4  >  z  I 
N1  VI  (5)  '  1 
MM  to)  '  1 

n  - 

11 
I  L 
\l 
A 

l.  AM 

I  UP 
I  At  A 
VO 
C.A 

OMI  li.AI 

omic.a: 

I’ 

0 

RIIO 


maximum  numher  of  mikles  to  he  re.ul  in  mote  that  the  program 
allows  the  nuinhei  ol  movies  reail  in  to  van  with  frequency  1 

uulieates  which  iiuantities  are  to  he  litteil  In  a  euhie  spline 
t\l  KM  Mi.  \  I  K  \Mi.  Kl  I  IIP.  IM  HIP.  \  ll  ami  PIIVOC  are 
moile  inputs  to  he  ilesenheil  subsequently  1 

then  a  euhie  spline  tit  is  applied  to  \1  KM  \(» 

then  a  euhie  spline  I  it  is  applied  to  \  I  K  AM  1 

then  a  euhie  spline  tit  is  applied  to  Kl  l  HP 

then  a  euhie  spline  lit  is  applied  to  l\l  I  IIP 

then  a  euhie  spline  lit  is  applied  to  PH \  Ol 
then  a  euhie  spline  tit  is  applied  to  \  I  I 

is  the  numher  of  integration  intervals  in  the  fiequency  range 

dill  1 

upper  frequency  ol  integration  in  kilohertz, 
lower  frequency  of  integration  in  kilohertz 
numher  of  frequencies 

an  array  of  four  elements  used  to  describe  the  source  function 
given  111  1  q.  t").  I  he  units  of  \  are  amperes 

an  arrav  of  tour  elements  used  to  describe  the  source  function 
given  in  1  q.  t  I  he  units  ot  ll  AM  t>  in  the  equation)  are  inverse 
seconds. 

characteristic  tune  in  seconds  associated  with  the  source  function 
given  m  I  q.  t 

characteristic  time  in  seconds  associated  with  the  source  function 
given  in  I  q.  1 7). 

characteristic  velovitv  in  m  see  associated  with  the  source  function 
given  in  1  q.  (  '). 

gain  for  the  receiver  function  given  in  l  q  ti'l 
angular  frequency  for  the  receiver  function  given  in  I  q.  (0) 
angular  frequency  for  the  receiver  function  given  in  1  q.  ltd. 
integer  variable  used  in  the  receiver  function  given  in  l  q  (<d. 
integer  variable  used  in  the  receiver  function  given  in  1  q.  tot. 
transmitter-receiver  distance  in  kilometres  used  in  mode  sum 


s 

S  =  1  for  positive  Fourier  transform 

S  =  -1  for  negative  Fourier  transform. 

INTPRT 

(lag  to  control  the  print  interval.  The  first  20  values  are  printed 
and  then  every  INTPRTth  one.  For  example  if  INTPRT  =  10  then 
the  1  st  20  values  are  printed  followed  by  the  30th.  40th.  50th  .  .  . 
every  1 0th  value  out  to  the  end. 

IP  LOT 

Hag  to  determine  whether  or  not  plots  are  drawn. 

If  IPLOT  =  0  no  plots  are  generated. 

It  IPLOT  =  1  six  plots  are  generated:  source  spectrum,  receiver 

spectrum,  channel  spectrum,  product  spectrum  (source*receiver* 

channel),  output  waveform,  and  the  input  current  pulse.  a 

TMIN 

an  array  of  three  elements  used  to  describe  the  starting  time  in 
seconds  for  the  input  current  pulse  plot. 

TINC 

an  array  of  three  elements  used  to  describe  the  time  increment  in  ’ 

seconds  lor  the  input  current  pulse  plot. 

NUMTS 

an  array  of  three  elements  used  to  describe  the  number  of  times 
that  are  plotted  on  the  input  current  pulse  plot. 

TAUMAX 

controls  the  latest  time  in  seconds  plotted  on  the  output  waveform 

curve. 

The  mode  data  or  third  part  of  the  input  follows  the  namenst  input.  The  eight  col¬ 
umns  of  mode  data  on  pages  9  and  1  0  are: 

NMF 

1 

number  of  modes  at  FRHQ( I).  1  =  1.2 . NF  (column  1  I. 

RLTHPt  M.I) 

the  real  part  of  the  complex  ground  eigenangle  for  mode  M  and 
frequency  1  in  degrees  (column  2). 

IMTHPt  M.I) 

the  imaginary  part  of  the  complex  ground  eigenangle  for  mode  M 
and  frequency  l  in  degrees  (column  3). 

XTRMAC.(M.I) 

magnitude  of  excitation  factor  for  mode  M  and  frequency  1  (col- 
until  4). 

XTRANGt  M.I) 

phase  (in  radians)  of  excitation  factor  for  mode  M  and  frequency  1 
(column  5). 

FRF.Q(I) 

frequencies  in  kilohertz  for  which  mode  data  is  input  (column  6). 

ATT(M.I) 

attenuation  rate  in  dB/1000  km  for  mode  M  and  frequency  I  (column  7). 

PHVOC(M.I) 

phase  velocity  over  free-space  velocity  for  mode  M  and  frequency  1 
(column  8). 

The  mode  data  input  shown  on  pages 9  and  10  is  fora  single  mode  case.  A  sample  input  of 
mode  data  for  a  multimode  case  is  shown  on  pages  1  1  and  1  2.  The  eight  columns  have  the 

same  meaning  as  above. 

The  ordering  is  such  that  all  modes  for  the  first  frequency  are  fol- 

lowed  by  all  modes  for  the  second  frequency,  etc.  It  should  he  mentioned  that  the  attenua¬ 
tion  and  phase  velocity  inputs  are  not  used  directly  in  the  calculations.  They  are  included  in 
the  input  so  that  they  may,  at  the  users  option,  he  spline  fit  for  the  purpose  of  explicitly 
exhibiting  their  frequency  dependence. 
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IV  PROGRAM  LAYOUT 


This  section  describes  the  basic  features  of  the  pulse  shape  program  listed  in  the 
appendix. 

Reading  and  printing  of  input  quantities  occurs  in  MAIN.  MAIN  calls  in  order  the 
following  subroutines. 

SUBROUTINE  I  UNSPL  <MI>,LF.XX,YY,B.(  ,1)1 

Inputs  to  IUNSPL.  are  a  mode  index  MO.  which  takes  on  values  I  through  NM.  and 
the  index.  LT.  for  the  quantity  which  is  to  he  approximated  as  a  function  of  frequency  by  a 
cubic  spline.  LI  can  take  on  integer  values  I  through  h.  I  UNSPL  calls  the  two  following 
subroutines. 


a.  SUBROUTINE  FUNCVT  (MD.XX.YY)  places 

XX(i)  =  FRFQ(I).  1=1.: . NT. 

YY(k>  =  XTRMAtitMD.K )  it  1  F  =  1  and  data  read  in  for  l=k. 

YY(K)  =  X  TR  AN(i<  MD.k )  it  I  T  =  :  and  data  read  in  for  l  =  k. 

YY(K  I  =  RFT1  IP( MO.K )  if  LT  =  3  ami  data  read  in  for  l  =  k. 

YY( k )  =  IM'TIIPtMD.k )  it  LT"  =  4  and  data  read  in  for  l=k. 

YYtkl  =  PHVOCtMO.k)  if  l.T  =  5  and  data  read  in  for  l=k. 

YY(  k  t  =  A  I"T( M O.k I  it"  l.T"  =  0  and  data  read  in  for  I  =  k. 

b.  SUBROU  TINT  SPI  INI  (XX.YY.B.C’.D.N)  determines  the  coefficients  B.C.  O. 
of  a  cubic  spline  interpolating  the  given  curve  <  X.X( It.  YY(I).  I  =  1 .  2.  .  .  .  N ).  If  XXtll.LF. 
X.Lh.XXt  1+ 1 )  and  1 1  =  X-\\(  I ).  then  the  interpolated  value  at  X  is  Ft  X)  =  YY(  I )  +  B(l)  * 

II  +  U(1)  *  11  **  2  +  0(1 1  *  11  **  3.  The  interpolated  value  is  evaluated  using  the  function 
SPF VAL( XVAL.X.Y, B.C. O.N. INTI").  In  particular  SPTVAL  evaluates  the  interpolating  cubic 
spline  for  the  data  ( X(  I).  Y(  I )),  I  =  I  ...Nat  XV  AT.  INTI"  is  an  estimate  of  the  interval 
where  XVAL  lies.  X(  IN  IT  1. 1  F.XV  AL.LF.Xt  IN  IT  +  1 1.  but  need  not  be  used .  Set  IN  IT  =  0  if 
there  is  no  estimate.  On  return,  INIT  will  contain  the  interval  number. 

The  replacements 

Y(’(  LF.MO.I)  =  YY(I) 

BC  (LF.M0.ll  =  B(l) 

CCtLF.MD.I)  =  U(0 
DO  LF.MO.I)  =  0(1) 

are  then  made  in  FUNSPI  and  control  returned  to  MAIN.  MAIN  then  calls  the  TRXMTR. 
RTC'VR  and  (Ml  ANT  1  subroutines  at  the  frequency  points  F  =  ( k-  1 )  FU-FL  2^  +  FL. 
k  =  I.  2  .  .  ..  2n  +  I. 

SUBROUTINI  TRXMTR(F)  calculates  the  spectrum  for  the  Williams  description  of 
the  mean  lightning  stroke  as  discussed  in  Section  II.  The  subroutine  can  be  easily  altered  to 
satisfy  the  users  needs.  Tor  example,  alternative  descriptions  of  the  lightning  discharge  are 


readily  accommodated.  In  the  example  in  Section  VI  of  this  report  illustrating  results  for  a 
spread-spectrum  system,  the  transmitter  function  given  by  Kq.  (20)  of  that  section  was  used. 

SUBROUTINE  RFCVR(F)  calculates  the  spectrum  for  the  receiver  function  discussed 
in  Section  II.  This  subroutine  is  also  easily  altered  to  satisfy  the  user’s  needs.  For  example, 
the  spread-spectrum  system  calculation  presented  as  an  example  in  Section  VI  of  this  report 
utilized  the  receiver  function  given  by  Hq.  (21 )  of  that  section. 

SUBROUTINE  CHANHL(F)  calculates  the  spectrum  for  the  elf/vlf  channel  described 
by  Eq.  (4)  of  Section  II.  Specifically  it  is  for  the  vertical  electric  field  at  the  ground  produced 
by  a  ground-based  vertical  electric  dipole.  Ground-based  horizontal  dipole  sources  can  be 
accommodated  using  the  replacement  indicated  by  Eq.  ( 12)  of  Section  II. 

The  real  part  of  the  product  spectrum,  Idl(F) R(F)ll(F,p),  which  occurs  in  Eq.  ( I )  of 
Section  II  is  stored  in  X(k)  and  the  imaginary  part  is  stored  in  Y(K). 

SUBROUTINE  NLOGN(N.X.Y.SIGNT.A.B)  calculates  (apart  from  end  point  effects) 
integrals  of  the  form  (S  =  SIGNT) 

f" 

exp  l  — j 27rS At |  I  (X(F)  +  j Y ( F ) )  exp(j27rSFr)dF 

J  A 


•b-a 


(X(E+A)  +  jY(f'+A))  expCi 2rrSFr)dF 


(14) 


by  the  fast  Fourier  transform  technique  of  Cooley  and  Tukey  (Ref.  4).  This  makes  use  of 
digital  evaluations  at  the  frequencies 


F(L)  =~~T7"  (B-A):  L  =  1.2  ..  ..  2N 
->N 


(15) 


and  the  method  yields  evaluations  for  the  times 


r<  K ) 


K  -  I 
B-A 


:  K  = 


-iN 


(lb) 


Real  and  imaginary  parts  of  the  integral  are  then  stored  in  X(K)  and  Y(K)  respectively. 
NL(X!N  also  has  built  into  it  the  Filon  weight  factors 


4(B- A)  2N 

s  1  * 

(K  -  I )-  (27U- 


2tt(K  -  I ) 
->N+1 


;  K  = 


1.  2.  5 


sN 


(17) 


which  for  K  =  1  is  simply  the  integral  size  <B-A)/2  .  If  the  integrand  of  Eq.  (14)  at  the 
points  (A.B)  is  not  negligible,  it  is  necessary  to  add  to  X(K)  and  Y(K)  the  following  end 
point  corrections 


B-A  [~  I  t 
2tt(  K  -  I)  j(SIGNT) 


sN 


2tt(  K  -  I) 


( I  -exp(-jS27r(k-  1)/ 


:N„][- 


U(  1  ,K)  +  U(2N  +  1  ,K) 
(IS) 


14 


where  the  U\  .ire  the  complete  mteni.mil  of  Hq.  1 14)  evaluated  at  the  frequencies  |see  Hq. 

1 1 5 ) |  I  I  and  I  2^  ♦  I  It  K  I  in  the  factors  multiplying  the  U’s  in  Hq.  ( IS)  the  multi¬ 
plying  factors  become  one  It  ill  the  interval  si/e,  (B- A)/2N+1 .  If  the  Filon  weight  factors 
were  omitted,  the  weight  laclors  given  by  hq.  (  17)  would  simply  be  replaced  by  (H-A)/2^ 
and  I  q  (IX)  by 

ii— -  f-l  t  l.k)  +  U|2N  +  I.Ktl  (Id) 

sN+l  L  J 

l  he  quantity  S  =  SKINT  takes  on  the  values  +  I  or  -  I  and  simply  allows  for  plus  or  minus 
transform  as  desired.  It  should  be  observed  that  although  the  region  of  significance  of  the 
integrand  of  I  q.  1 14)  may  be  quite  small,  the  integration  limits  A.B  may  of  necessity  be 
quite  large  in  order  to  achieve  a  desired  time  resolution  |see  Hq.  ( l6)|.  N  must  he  chosen  to 
give  small  enough  step  si/es  in  the  region  where  the  integrand  is  significant.  Specifically,  step 
si/es  must  be  small  compared  with  distances  (in  frequency  units)  over  which  the  integrand 
(exclusive  of  the  exponential  factor  when  Hilon  weight  factors  are  used)  changes  appreciably. 
Also,  it  should  be  mentioned  that  the  program  can  be  easily  altered  to  accommodate  other 
integration  routines  should  the  need  arise  in  a  particular  application. 


V  DHSURIPTION  OF  OUTPUT 

The  sample  output  shown  on  pages  1  7  through  24  begins  with  a  listing  ol  the  three 
plot  identification  lines  followed  by  the  namelist  output.  The  mode  data  come  next.  For 
each  frequency  (given  in  increasing  order)  the  number  of  modes,  the  real  and  imaginary  parts 
of  the  eigenangle,  the  magnitude  and  phase  of  the  excitation  (actor,  the  attenuation  rate  and  the 
phase  velocity  normalized  to  tree-space  velocity  are  listed.  I  hough  the  sample  output  is  lor 
a  single-mode  case,  the  program,  as  mentioned,  is  equally  suited  for  multimode  studies. 

The  principal  output  of  the  pulse  shape  program  begins  on  page  Id.  The  transmitter, 
receiver,  channel,  and  product  (XMTR*RC'VR*CHNL)  spectra  are  given  as  a  function  of  fre¬ 
quency.  Not  all  2044  (i.e.,  2^  +  I  with  N  =  I  I )  lines  are  listed.  The  printout  is  controlled 
by  the  namelist  variable  INTPRT.  The  first  20  values  of  the  spectra  are  always  printed.  These 
are  followed  by  every  20th  value  of  the  spectra  because  in  this  instance  INTPRT  =  20. 

Following  the  spectra  output  comes  output  (page  22)  pertaining  to  the  time  signature 
of  the  output  waveform,  (i(p.t)  given  by  Hq.  1 1 ). 


The  first  column  is  the  time  in  seconds,  the  second  and  third  columns  are  the  real  and 
imaginary  parts,  respectively,  of  the  integral 


i  dl(cu)r(  w)h(co.p)  expljcot  )dw 


The  last  column  is  the  waveform,  ( «< p, t ),  given  by  Hq.  ( I )  in  volts/m.  CKp.t)  at  time  equal  to 
zero  should  be  zero.  The  departure  from  zero  at  time  equal  to  zero  is  believed  to  he  associated 


* 


wi'li  truncation  effects  ami  or  discontinuities  in  the  third  and  higher  derivatives  of  the  inter¬ 
polated  mode  data.  I  he  program  also  generates  six  plots.  These  are: 

1 )  transmitter  spectrum  vx  freij 

2)  receiver  spectrum  vx  Iren 
.))  channel  spectrum  vs  Iren 

4i  product  (transmitter*receivcr*channcl)  spectrum  vs  fren 
5)  output  waveform  (i  vs  lien 
0)  input  current  pulse  vs  fren 

l  he  four  spectra  plots  are  only  plotted  between  the  first  and  last  frequency  inputs 
<1  HI  0<H-l  RKXNI  )l 

The  output  waveform  is  plotted  out  to  TAUMAX  (a  namelist  input  variable), 
l  he  plots  are  shown  in  Figs.  I  through  (>. 
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-  .253' C-007 
- . 75296-008 

.20401-007 


62633*000  .20096-007 
63300*000  .27444-097 
63907*000  .117.1',!  007 
64633  *  000  -.46412-008 
66300*000  -.13730-007 
66907*000  -.66636-008 
66633*000  .92403-009 
67300*000  .26990-009 
67967*000  .13600-011 


.  95583-009  .  4  t  7 ' M  - 0 0 7 
.16196-007  .44000-007 
.24031-007  .23497-007 
.20346-00/  -.92023-000 
.15976-007  -.27460-007 
.27011-000  -.13107-007 
.69037-000  .10401-000 
.53300-000  .53996-009 
.10206-008  .27216-011 


INPUT  CURRENT  PULSE  (AMPS) 


1 

- 1 - » 

24.0  32.0 

TAU  (MILLISECONDS) 

TRANSMITTER: 

WILLIAMS  SOURCE 

RECEIVER: 

P=2,  Q=2,  F1  =  10  HZ,  F2=2500  HZ 

CHANNEL: 

SATELLITE  NIGHT  A=254,  C=47,  RHO=3700  KM 

Figure  5. 

Output  waveform. 

8 

O 

CM  -1 
"1 

8  l 

WILLIAMS  SOURCE 

40  60  80  100  0  200  400  600  800  1000  0 


TIME  (MICROSECONOS) 


TIME  (MICROSECONDS) 


2  4  6  8 

TIME  (MILLISECONDS) 


Figure  6.  Input  current  pulse. 


VI  AN  ADDITIONAL  APPLICATION 


As  .in  .ulilition.il  example  ol  the  problem  type  to  wlueli  the  present  program  may  he 
applied  we  consider  a  ease  study  similar  to  that  examined  In  Rotlimuller  (Ref.  5)  in  his 
teport  on  the  etleet  ol  the  propagation  channel  on  spread-spectrum  communication  systems. 
Rothmuller  investigated  the  etleet  that  propagation  at  very  low  frequencies  (vIO  through  the 
earth-ionosphere  waveguide  has  on  one  type  of  communication  system.  I  he  system  studied 
v»as  characterized  by  a  dilterenti.il  phase-encoded  signal  waveform  composed  by  frequency 
shift  keying  ( I  SK )  a  carrier  with  a  binary  pseudo-random  or  pseudo-noise  (PN)  sequence  of 
puls*,  s  or  *  hips.  I  lie  LSk  mod  illation  index  is  0.5 ,  which  is  designated  as  minimum  shift  key¬ 
ing  (MSK).  Lor  more  detail  concerning  the  basic  waveform  and  terminology  the  interested 
reader  is  referred  to  Kothnuiller's  report.  Mere  we  note  only  that  the  PN  sequence  lias  a 
power  spectrum  given  by 


where 

l'0  is  the  carrier  frequency 
b'c  is  the  chip  frequency. 


I  he  communication  system  to  be  evaluated  is  assumed  to  consist  in  part  of  a  receiver 
followed  by  a  demodulator  matched  to  the  undistorted  transmitter  signal.  The  receiver 
response  modeled  by 


.+jF_i!o  3 

F, 


F  +  F  I  3  :  F,  =  I  kHz 

I  -  j - - 


will  be  assumed  in  the  subsequent  calculations. 


Figure  7  shows  a  vlf  waveguide  signal  plot  as  a  function  of  frequency  for  a  daytime 
Hawaii  to  southern  California  propagation  direction  and  for  a  path  length  of  2 28 2  km. 
Observe  the  deep  null  at  23.5  kll/.  One  possible  measure  of  relative  performance  (RP)  of  a 
spread  spectrum  system  operating  at  a  carrier  frequency  Ft,  =  23.5  kll/.  over  the  same  system 
operating  at  the  single  frequency  F0  is 
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50  I— 
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.* 
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HAWAII  SOUTHERN  CALIFORNIA  DAYTIME  PROPAGATION  PATH 
P  -  2282  km 


23  24 

FREQUENCY. kH; 


Figure  7.  Daytime  mode  sum  as  a  function  of  frequency 


where  the  subscript  “max"  signifies  the  maximum  value  of  the  squared  denominator  I  he 
Fourier  integrals  in  Hq.  (22)  can  be  evaluated  using  the  present  program  and  results  are 
shown  in  Fig.  8  for  chip  frequencies  of  100  sec-'  and  1000  sec-' .  It  will  be  seen  that  two 
correlation  peaks  occur.  This  phenomenon  has  been  discussed  by  Rotlunuller.  A  relative 
performance  of  20  dB  could  be  expected  for  the  case  of  F\-  =  1000  sec-'  and  about  4.5  dB 
for  the  system  operating  at  Fc  =  100  sec-' .  Of  course,  if  the  system  were  operating  at  a 
central  frequency,  F0,  corresponding  to  a  maximum  in  the  mode  sum,  the  relative  perform¬ 
ance  would  be  degraded.  Generally,  though,  the  gain  in  performance  in  the  neighborhood 
of  nulls  would  outweigh  the  loss  of  performance  in  the  neighborhood  of  maxima. 
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IMPLICIT  REAL -81 A-H.O-Z ) 

COMPLEX. 16  IEMPf.TEMPL.GB 
COMPLEX. 16  XMT R , RCVR , CHNL . FOF  TAU 
COMPLEX* 16  1M/10.0D0.1 .ODO)/ 

REAL. 8  IT 
RE  A  L  *  8  IMTHP 

REAL* A  PLOTX| 3000) , P  LOT Y 1 (3000) , P LOT Y2 ( 3000 ) , P LOT Y3 ( 3000 ) 

REAL'A  PLOTY  4  I  3000 ) 

INTEGER  P,Q 
COMMON  ONE/XMTR 
COMMON/ T  WO/ RCVR 
COMMON/TMREE  CHNL 

COMMON,  FOUR  A ( 4 ) , GAM ( 4 ) , T AUP , T AUV , VO 
COMMON  FIVE,  GA .OMEGA  1  , 0MEGA2 , P , Q 
COMMON, S I  X 'MODE (15.50) 

COMMON  SEVEN  F  HE  0 ( 50 )  . XT  RMAGt  15,50),XTflANG(15.50),RETHP(15.50). 

$  IMTHP l 15,60) .ATT ( 15.50) , PHVOC ( 1 5 ,50 ) .NFVF.NF.NM 
COMMON  NINE  WHO 
COMMON  E  LEVI N  K K ( 15) 

DIMENSION  TMIN( 3) , T INC( 3) ,NUMTS(3) 

OIMC.NS  ION  XX  (  50  )  ,  YY(  50  )  .  B(  50)  ,  C(  50),  D<  50) 

DIMENSION  NEVF ( 6) , X(2049 ) , Y( 2049) 

DIMENSION  LABEL! (10), LABE  LR ( 10), LABELCI 10 ) 

NAMEL I  ST /DAI UM/A , GAM, TAUP, TAUV . VO.GA .OMEGA  1  , OMEG A2 , P , 0 , NM , NE V F . N . 
$  FU.FL.NF .RHO, S, INTPRT,  IP  LOT , TAUMAX 
$  , TMIN, T INC .NUMTS 

DATA  PI/3. M I59265358979324D0/ 

READ( 5,901 )  LABELT 
PRINT  902. LABElT 
READ( 5.901)  LABELR 
PRINT  903, LABELR 
RE AD( 5 . 90 1 )  LABELC 
PRINT  904 , LA  BE  LC 
READ! 5, DATUM) 

WRITE ( 6. DATUM) 

TWOPI  =  PI-2. ODO 
DTR  =  PI/1 .8002 
DO  25  K  F  a  1  ,  N  F 
DO  25  M  =  1  ,  NM 
M0DE(M,KF)=0 
PRINT  920 
DO  92  KFsl  ,  NF 
M  =  0 
NMS  -  0 
M  =M.  1 

READ! 5 . 1 1  I  NMF . RE  T HP ( M , N F ) ,  IM THP( M , K F ) . XTRMAG (M.KF) , XT  RANG1M , KF ) , 
$  FREOI KF ) , AT T (M.KF ) , PHVOCl M, KF ) 

MR  I TE ( 6. 1 2 '  NM' . F  R EQ ( KF )  , RE TH P ( M . KF )  , IMTHP (M.KF)  . XTRMAG (M.KF ) . 

S  X  TRANG(  M.KF  )  .ATT  (F.1.KF)  .PHVOC  (M.KF) 

IF  (NMF  .NE.  0)  NMS=NMF 
I F ( NMS  .  EO.  0)  NMS'NM 
RETHP(M.KF)3RETHP(M,KF)*DTR 
IMTHP (M,KF);IMT HP (M,KF)*DTR 


MODI  iM.KFliVI 
I  F  |  M  .LI.  N’.'s  )  GOT 09 1 
92  F RIOt kF  |  .FRt OIKF )• 1000 .00 
F  L  -  F  L •  1000.00 
FU  FLJ  •  1  000 . 1)0 
00  7(,  M  «  1  ,  NM 
KMOI'l  •  0 
KM"L'I  1  1 

00  ! 1  KF  •  1  ,  NF 

71  IF (Moot (M.KI  I  .NE.  0)  KMODE  *  KVOOE ♦ 1 

00  -1  Kl  •  I  .  Nl 

74  |  F  l  MODE  (  M  .  K*  I  .  NE  .  01  MODE  l  M  ,  KF  )  -  KMOOE 

DO  KF.1.NF 

IF iMODElM.KF  )  . EJ .  01  KMOOE 1  *  KMOOE 1 ♦ 1 
KK  l  M '  *  KMODt  1 

75  IMVOOEIM.KE  I  -NE.  0)  GO  TO  70 

76  CONItNUE 

DO  65  MD-1.NM 
LF  *0 

L  0 

50  l’l)l 
NFV I  •  L 

51  IF  (NEVFUI  .EO.  1)  GO  TO  59 

L  F  *  t  f  ♦  1 

CALI  FL'NSPL  IMD,  LE  ,  XX, YY  .  B.C,  0) 

59  1 F  i L  . LT.  6)  GO  TO  50 

65  COM  1NUE 
NN  =  0 
NFV.'.  .  N 
NF'1  NP*1 

Dl  T  (  ru-Fl  1  NP 

P  f;  ini  “15 

NUMPTS  >  0 

DO  10  K  1  .  NP  1 

F  M  K  II • ( FU- F L 1  NP+F  L 

CALL  1 RxMTR  i  F  1 

CALI  WE  CVM(  F  1 

CAL  L  L'FIANC L  (FI 

X  I  K  I  XMT  H  •  Hi  V  R  *  CUN  L 

Y(K'  -  IM«  V.'l  K*RCVR*CfiNL 

1 F  L  f  .LI.  FtUO(')  .OH.  F  .  GT  .  FREOlNFI)  GO  TO  6« 

NL'MF'  1  5  =  NUMP  1  •>*  1 
PLOTX  l  NUMPT'.  I  »  F 

PlOTy 1 INUMPI 91  ■  20.0  *  DL0G1 0(CDAD9(XMTR1 ) 

PLOTy  2  (  NLIMP 1  9  1  .  20.0*DLOG10|  CDABSiRCVRI  1 
PLOT*  JlNlIMPt  9  I  •  20 . 0  •  0L0G1  0  (  COABS  L  CHNL  1  1 
PLOt  >  4( NUMP1 s 1  20.0  *  D10G1 0 ( D90R  T (X(M**2  +  Y(K)*»21) 

66  | F  ( K  .Gl  .  20  .AND.  M00( K .  I  NT  PR T  )  . NE .  01  GO  TO  10 
PH  ’.  Nt  900  .  F  .  \MTH  ,  RCVR  ,CHNl .  XL  K  1  ,  V  (K)  .  K 

10  CONTINUE 

l F (  t  P  LOT  .EG.  01  GO  TO  250 
CALL  1NK(  '  Pt  N  2  »  BLACKS'  1 
CALL  BGNPLL1) 

CALL  CPNPLT 

CALL  XLABELt ' FRE01M21 1 ,61 

CALL  V  LABEL!  ' 20*  LOG  1 0 1 ABU l XMT  R ) ) ' ,19) 

CALL  XLGPLT(PL0TX.PL01Y1 .NUMPTS) 


1 


call  MESSAGl  ' TRANSMITTER  SPEC T RUM : • , 2 1 , 2 . 0 , 8 . 5 ) 

CALL  MESSAGl  LABELT  ,40,2.0,8.3) 

CALL  EMDPL(O) 

CALL  BONPLUI 
CALL  OPNPLT 

CALL  XL  ABEL!  ' FREQ|HZ )  '  ,8  I 

CALL  V  L  ABE  L (  2O*LOG1O(A0S(RCVR))’,19) 

CALL  XLGPLT  l  PLOT  X . PLOT Y2 .NUMPTS) 

call  MESSAGl 'RECEIVER  SPEC T RUM : 1 , 1 8 , 1 . 0 , 9 . 6 ) 

CALL  MESSAGl LABELR, 40, 1 .0,9.4) 

CALL  ENDPLtO) 

CALL  BGNPLIJI 
CALL  OPNPLT 

CALL  X  LABE  LI  ' FREQ(HZ) '  ,8) 

CALL  Y  LABE  LI  ' 20- L0G1 01 ABS( CHNL  )  )  *  .19) 

CALL  XLC.PLT  |  P LOT  X, PLOT Y 3. NUMPTS) 

call  MESSAGl  1  CHANNEL  SPEC T RUM :  '  ,  T 7 , 1  . 0 , 8 . 4 ) 

CALL  MESSAGl LABELC, 40, T . 0 ,8.2 ) 

CALL  ENDPLlOl 
CALL  BGNPL ( 4 ) 

CALL  OPNPLT 

CALL  XLABELl ' FREQ(HZ) ' .8) 

CALL  YlABEU  1 20*L0G10| ABS( XMT R • RCVR* CHN L ) ) 1 ,29) 

CALL  XLGPLTIPL0TX.PL0TY4, NUMPTS ) 

CALL  MESSAGl ' PRODUCT  SPEC T RUM ' . 1 6 , 0 . 5 , 1 0 . 0 ) 

CALL  MESSAGl 'TRANSMITTER:', 12, 0.5, 9. 6) 

CALL  MESSAGl  LABELT  ,40,2.0,9.6) 

CALL  MESSAGl  'RECEIVER:', 9, 0.5, 9. 4) 

CALL  MESSAGl LABE LR , 40. 1 . 7 ,9 • 4 ) 

CALL  MC  SSAC.I  '  CHANNEL  :  '  ,8, 0.5, 9. 2) 

CALL  MESSAGl LABELC, 40. 1 .6,9.2) 

CALL  E  NDP  L| 0  I 
250  CONTINUE 

TEMPF  «  X (  1  )♦ I M» Y( 1  ) 

TEMPL  =  X | N  P 1  ) ♦ I M • Y ( N  P 1  ) 

CALL  N  LOGN |N,X,Y,S,FL,FU) 

PRINT  905 
NUMPTS=0 
DO  20  K  =  1  ,  Np 
TAU=|  K— 1  )  I F  U-F  L ) 

OM  =  T  WOP  1  *  T  AU 
I F i K  . NE .  11  GO  TO  2 
GB  =  OLT  2.000 
GO  TO  3 

2  GB  DLT  |S* I M • OM  ) ♦ ( 1 . ODO-CDE X P ( - IM» S *0M  * DL T ) ) / | OM*OM ) 
GB  =  GB  DLT 

3  FOFTAU  -  CDE.  XP  l  I V*  S  *0M  *  F  L  )  *  |  X  |  K  )  +  IM*  Y  |  K  )- T  EMPF  .  GP  + 

$  T EMP L • GB *CDE XP ( IM»OM  ♦S*(Fu-Fl))) 

RFOFT  =  2. 0D0*DREAL( FOFTAU) 

I F |  T  A  U  .GT.  TAUMAX)  GO  TO  24 

NUMPTS  *  NUMPTS+1 

plot x i numpts *  =  tau*i.0E3 

PLOTy 1 (NUMPTS)  *  RFOFT 
24  CONTINUE 

IF|K  , GE .  20  .AND.  MODI K , INTPRT )  ,NE.  0)  GO  TO  20 
PRINT  910. TAU, FOFTAU, RFOFT 
20  CONTINUE 


t 


0)  GO  TO  260 


IF  l  l  P  LOT  .fO. 

GALL  BGNPL<5) 

CALL  OPNPLT 
CALL  XL  ABF  L(  '  T AU( MI LL I  SECONDS )  '  , 1 7 ) 

CALL  YLABELI  'OUTPUT  M A V E F ORM-G ( V ,  M ) *  . 22 ) 

CALL  L INPl T( PLOTX , PLOTY 1 . NUMPTS ) 

CALL  MESSAGE  1  OUTPUT  WAVEFORM'  , 15,0.0.8.8) 

CALL  MESSAGl  TRANSMITTER:  '  , 12,0.0,8.4) 

CALL  MESSAGl  LABELT  ,40,1.5,8.4) 

CALL  MESSAGl  RECEIVER: '.9. 0.0, 8. 2) 

CALL  MESSAGl LABELR. 40,1 .2.8.2) 

call  vessagi ‘channel: 1 .b.o.o.b.o) 

CALL  MESSAGl LABELC, 40 . 1 . 1 ,8.0 ) 

CALL  E  NOP  HOI 
CAlL  BGNPC16I 
CALL  INTAXS 
T  =  T  M  I  N  t  1  t 
DO  45  J*  1  , NUMTSE  1  ) 

'.  T  3  0 . 000 
DO  40  A  =  1  .  4 

IT  IT«A(HI*DEXP(-T *GAM ( K ) ) 

40  CONTINUE 
PLOTXIJ)  =  1  • 1  . OE  6 
P  LO  T  >  IIJI  =  IT 

T  s  T ♦ T I NC l 1 ) 

45  CONTINUE 

CALL  PHYSORi C .8 .2 .0  I 

CALL  T  I  T  l E  t  '  -1  ,'  T I ME ( MICROSECONDS) 18. 

$  'INPUT  CURRENT  PULSE | AMPS ) ’ . 25 . 3 . 0 . 5 . 0 ) 

CAlL  GRAF(0.  .20. ,100.  .0.  ,3000.  ,  15000.  ) 

CALL  CURVE! PLOTX, PL0TY1 , NUMTSl 1 ) , 0 ) 

CALL  endgreoi 
T  =  T  M  I  N  (  2  ) 

DO  46  J=1 .NUMTSI2) 

IT  3  0.000 
DO  4 1  K  -  1  ,  4 

IT  3  1  TxA|KI  -DEXPI-T^GAMIK)  ) 

41  CONTINUE 

P LO  t  X  I U  )  3  T -  1  . 0E6 
P  L  0 1  Y 1 ( J I  3  IT 
T  3  T ♦ T I NC ( 2 ) 

46  CONTINUE 

CALL  PHYSORi 4.25,2.0) 

CALL  TITLE!'  ' ,-T , ' TIME(MICROSECONDS) 1 . 18. ’  ’,1,3. 0,5.0) 
CALL  GRAF | 0 . .200. ,1000. ,200. .400. ,22 00.) 

CALL  CURVE!  PLOT  X.  PL0TY1  ,  NC'MTSI  2)  .  0  ) 
call  MESSAGl LABELT ,40, 1 .0,5.6) 

CALL  ENDGR(0l 
T  =  T  M  I  N  i  3  ) 

DO  47  J3 1 .NUMTSI 3 ) 

IT  3  O.ODO 
DO  42  K= 1 ,4 

IT  3  I T+AiK) • DEXP( -T*GAM(K) ) 

42  CONTINUE 

PLOTX I J )  3  T • 1 . 0E3 
PLOTYUJ)  .  IT 
T  •  T+TINC(3) 


*4. 


47  CONTINUE 

CALL  PhySORI 7 . 70 , 2 . 0 ) 

CALL  TITLE!  ' .-1 , ' T IME( MILLISECONDS) ' , 18, '  ‘.1,3. 0,5.0) 

CALL  GRAF ( 0 . .2. , 10. . 100. , 100. ,600. ) 

CALL  CURVE l PLOIX , PLOTY 1 , NUMTSl 3 ) . 0 ) 

CALL  ENOGR(O) 

CALL  ENDPLlO) 

260  CONTINUE 

1 1  FORMA  M  15,2*  10.5.D15.6.4F10.5) 

12  FORMAT!'  '  .  15. 3F10.5.D15.6, 3F10. 5) 

900  FORMAT)  ‘  '.1CX.9D12.5.16) 

901  FORMA  T (  1 0A4 I 

902  FORMA  T (  '  1  '  , 1 0A4  ) 

903  FORMA  T  (  '  ‘  .  1 0A4  | 

904  FORMA  T l  *  '  . 1 0A4 .  ,  | 

905  FORMA  T('1'.46A,'TIME',7X,'RE(FFT)’,5X,'IM(FFT)',4X,'G(  RHO  , T ) -V/M '  . 
S  ,45X ,  '  ! St  vONOS )  '  I 

910  FORMAT! •  ' .42X.4012.5) 

915  FORMAT (' 1 ' , 1 JX ,' FRE0IH2 )' ,4X XMTR  R',6X.'XMTR  I',6X.'RCVR  R'.6X, 

$  1 RCVR  I  1  . 6  X ,  1 CHN  L  R '  , 6X ,  ' CHNL  I '  , 8 X ,  ' XMTR * RCVR • CHN L ' . 

$  1 0  X ,  1 K '  , / ,  '  ,98X,  'REAL'  ,8X,  '  1  MAG '  ) 

920  FORMAT !  '1  '  ,29X  ,  ' NMF  *  ,4X,  'FREO'  ,3X.  '  THETAR '  ,4X,  ’ THETAI ‘  ,7X, 

S  'XI R MAG '  ,7X,  'XTRANG'  .7X,  'ATT  1  ,6X, 1 PHVOC '  ,  /  . 36X ,  ' KHZ '  , 4X , 

$  'DEGREES' ,3X, 'DEGREES' , 19X, ' RADIANS' , 6X , ' DQ ' ) 

999  RETURN 
END 


SUHROUT  INE  MINSPL  (  MD  ,  LF  ,  XX  ,  YY  ,  B  ,C  ,  D) 

1MPL I C I T  REAL • BiA-H.O-/ ! 

Rl  Al  •  it  I  MT HP 

COMMON  SIX  MODE  (IS,  SO  ) 

COMMON,  SEVEN  I REQ<  50) , XTRMAOl  15,50), XT  RANG  I  1 5 . 50 ) . RE  T HP (  1 5 . 50 ) . 
S  IM1HP( 15,50 ), AT  T (  15,50 ) , PHVOCl 1 5 , 50 ) . NF V E . NF . NM 

COMMON  EIGHT, l M 

COMMON  TEN/1 Cl  6, 15 ,50) . BC ( 6 . 1 5,50  )  ,CC(6. 1 5 . 50 ) . DC ( 6 , 1 5 . 50 ) 

DIM!  NS  ION  XX| SO ) 

DIME  NS  ION  YM  50  )  . 

$  B|50  )  ,C(50  I  , D|50 ) 

CAl l  FUNCVF ( MO . XX . YY ) 

CALI  SPLINE! XX, YY.B.C.D, LM) 

DO  4 Li  1  =  1  .  NE 

YC ( LF , MD , I )  YY ( I ) 

DC ( LF ,MD, 1)  H( I ) 

CC (LF, MD , I ) - C ( I ) 

DC | LF ,MD, I )=D( 1 ) 

CON! 1 NUE 
RE  TURN 
END 


SU6W0UT  INI  I  i  NCU  |  Ml).  XX  ,  YY  ) 

IMP l  l Cl T  Wf  A l  MM  A-M.O-Z > 

W I  A  l  • 0  IMTMP 

coWi’N  s  i  x  ./onr  Mb.  ‘.'0 i 

COMMON  Si  VI  N  »  Rl  Ol bO) . XT  RMAGl  1  b  .  50  )  ,  XT  RANG  l  \  b  .  bO  )  ,  RETHP(  1  b  .bO  )  . 
$  IMTMPl Ib.SO) . ATT l 15.50) . PHVOC( 1 b . bO ) . NF V F . NF . NM 
COW  'N  I  K»MT  IM 
COY  -IN  IllVlN  nKi  1  b  ) 

01  Ml  NS  ION  XXlbO)  .YY(bO) 

GO  ID  ( J0.40.b0.60.70.80) .NFVF 
30  00  Jb  1  -  1  .  N» 

l  FIMODE l MO.  I  )  .10.  0)  GO  TO  3b 
JJ  -  l  KK  l  Ml)  1  ♦  1 

>  V i JJ  )  -  XTKMAGlMO ,  l ) 

X  X  l  JJ  )  *  »  Wl  01  1  ) 

LM  MODClMD.l) 

3b  CONTINUE 

GO  TO  “9 

40  DO  4b  l  1 .NT 

IFlMUDUMO.n  .10.  0)  GO  TO  4b 
J  J  I  -  KK  (  Ml)  »  ♦  1 

Y  > ( J  J  )  *  XTRANGlMD.  I  ) 

X  X  l  Jv'  T  T  •:  Oil) 

IM  MODE  l  Ml). -1  ) 

4b  CON!  I  Nut 

GO  TO  99 

bO  DO  bb  l  -  1  ,N» 

UlMODIiMD.l)  .10.  0)  GO  TO  bb 
0  J  1  -KM  MO  )  ♦  1 

>  Y (  v I v I  )  WE  1 H P l MD ,  1  ) 

XX  i  JJ  )  *  FRI  Ol  I  ) 

LM  MODE l MO, l ) 
bb  CONTINUE 

GO  Til  99 

60  DO  lib  l  1  .  N* 

l M MODE l MO, 1)  .EO.  0)  GO  TO  6b 
JJ  l - KK ( MO T ♦ 1 

Y  Y l JJ  )  -  l Ml  HP l MD ,  l  ) 

XX  t  J  J  )  f  Rt  Oil) 

IM  MODE l MO, 1) 

6b  CON  1  l NUl 

GO  TO  "9 

70  DO  7b  1-1  .  N» 

If  l  MODI  l  MD ,  l  )  .EO.  0)  GO  TO  7b 
JJ  I-KMM0W1 

Y  Y l J J  )  «  PMV OC l MD,  l ) 

XX(  JJ  )  -  FRI  Ol  1  ) 

LM  MOD El  MO, I) 

75  CONTINUE 

GO  TO  99 

80  DO  Ob  l » 1 ,N» 

1 F l MODE IMD, l )  .10.  0)  GO  TO  6b 
JJ  •  l -KK l MO ) ♦ 1 
YYi  JJ  )  ■  ATT IMD, I ) 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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SUBROUTINE  SPLINE  (X.  Y,  B,  C,  D,  N) 
IMPLICIT  REAL*8(A-H.O-Z) 

DIMENSION  X(t),  Y(t),  B('),  C(1),  D(  1  ) 


SPLINE  DETERMINES  THE  COEFFICIENTS  B,  C,  D. 

OF  A  CUBIC  SPt 1NE  INTERPOLATING  THE  GIVEN 

CURVE  ( XI  I  )  ,Y(  I  ) ) ,  1  =  1 . N.  IF 

X ( I ) . LE . XX . LE . X ( I ♦ 1 )  AND  H  =  XX  -  X(I), 

THEN  THE  INTERPOLATED  VALUE  AT  XX  IS 
F  l  X  X  )  =  Y  (  I  )  4  B I  I  I  *  H  +  C(  I  ) *H* • 2  +  D(  I  ) +H++3 . 
THE  INTERPOLATED  VALUE  CAN  BE  EVALUATED 
WITH  THE  FUNC I  I  ON  SP  EVAl. 

B.C.D.  MUST  HAVE  LENGTH  AT  LEAST  N. 

IF  (N.GT.2)  GO  TO  050 
Cl t I  =  0.0 
0(11  =  0.0 

B I  1  >  =  l  Y  I  2  )  -  Y  (  1  )  )  /  I  X  (  2  )  -  X  (  1  )  ) 

RE  TURN 

050  NN  *  N  -  1 
TB  =  0 

DO  100  I  =  1 , NN 
IF  IXIIH).LE.X(I))  GO  TO  800 
D(  I  )  =  X(  1  +  1  )  -  X  (  I  ) 

TA  =  ( Y( 1+1 )  -  Y ( I ) )  /  DII) 

Cl  I)  =  T A  -  TB 
T  B  =  TA 
TOO  CONTINUE 
Cl  II  =  0 
C  I  N  I  =  0 


T  A  =  0 
TB  =  0 

DO  200  I  =  2 , NN 

cm  *  cm  -  ta  ♦  cii-1) 

Bin  =  2.0  ♦  (Dm  +  D(I'I))  -  TA  *  TB 
TB  =  Dll) 

TA  =  TB  /  Bl  I  ) 

200  CONTINUE 


C I NN )  =  C ( NN )  /  B ( NN ) 

IF  (NN. LT .3)  GO  TO  350 
DO  300  I  =  3.NN 
J  =  NN  ♦  2  -  I 

300  CIJ1  =  (CUT  -  DIJ)  •  CIJ+11)  /  B(J) 
350  DO  -TOO  I  =  1  ,  NN 

B(I)  =  I  Y  (  I  ♦  1  )  -  Y  (  I  )  )  /  DID 

t  -  icui  +  cm  +  cmm  •  o m 
om  =  icu+n  -  cm )  /  om 
cm  *  3.o  ♦  cm 
400  CONTINUE 
RETURN 

BOO  PRINT  900 

PRINT  901  ,  I  , X ( I )  , X(  1  +  1  ) 


RETURN 

901  FORMAT  (IX.'  1  .',15,'  X(I)  *,.1PE12.5.'  X ( I +1 )  *'.1PE12.5  /) 

900  FORMAT  ('  ERPOR  IN  SPLINE',/, 

$  '  X-COOROINATE  VALUES  ARE  NOT  IN  INCREASING  ORDER') 

END 
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SUBROUTINE  TRXMTR(F) 

IMPLICIT  REAl.8( A-H.O-2) 

COMP LE X* 1 6  lM/tO.ODO. 1 .ODO)/ 

COMPLEX* 16  XMTR 
COMMON  ONE /XV T  R 

COMMON/ FOUR  A(4  ) ,GAM(4) , TAUP, T AUV , VO 
OATA  Pl/3. 141S92653B897932400/ 

TWOPI  *  P I *2 . ODO 
XMTR* (0. ODO, 0. ODO) 

OMEGA  =  T  WOP  I  *  F 
00  30  1  =  1,4 

XMTR=  XMTR+Al  I )/( ( GAM( I ) ♦ I M*OMEGA ) *  *2 ) •( 1 . DO-CDEXP 
$  (-TAUPM IM-OMEGA+GAM( I ) ))/(1 ,000+TAUV«( IM*OMEGA*GAM( 1)))) 
XMTR=XMTR-VO 
RETURN 
END 


SUBROUT INE  HtCVRl F ) 

IMP l IC1T  REAI •8IA-H.O-Z) 

COMP  L E  K* 16  1M  (O.ODO. t .ODO)/ 

C0MPLE\*16  RCVR 
INTIutd  P,0 
COMMON  TWO  RCVR 

COMMON  FIVE  GA .OMEGA t , OMEGA? . P , Q 
DATA  PI  3 . 1-11  I'92ob3b0lJ?9324D0  / 

TwOPI  *  PI*?. 000 
OMEGA *TWOPI*F 

RCVR'GA* { 1M*  OMEGA/OMEGA  1  )»*P/(  (  1 . 00+ IM*OMEGA/ 
(OMEGA  1  )**P.(  1 . 00+ 1M+0MEGA/ OMEGA? )**Q) 

RETURN 

END 
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SUBROUTINE  l MANE l ( F  ) 

IMPLICIT  REAL *8l A-H.O-Z) 

COMPi  E  X • 1 6  l ONST 
COM*' LEX*  16  1M  tO.ODO,  1  .000)/ 

COMPLEX* 16  MSUM.CPXSIN, EXC.CHNL 
REAL *8  I  M T  Hi' 

COMMON  THREl  CHNL 
COMMON  Six  MODE! 15.50) 

COMMON  SEVEN  1  REO(SO)  ,XTR.VAG(  15,50) ,  X T RANG (  15,50  )  ,RETHP(  15,50)  . 
$  1MTHPI  15,60  )  .  ATT  (15.50  )  ,  PHVOCl  1 5 ,50 ) . NF V F , NF , NM 
COMMON  NINE  RHO 

COMMON  TEN  u (6. 15,50) . BC 1 6. 1 5. 50) .CC I  6. 15,50 ) .DC ( 6. 1 5 , 50  ) 
COMMON  E  LEVIN  K  K.  I  15) 

DIMENSION  XX  I 50 ) . YY I  50 ) , 8(50 ) ,C( 50 ) , 0(50 ) , E(6) 

DATA  VL  I  TE  .  ‘19762505  ' 

DATA  PI  3.1-11 69205358S79324D0/ 

DAI  A  ERA?  i-  .  37103/ 

TmUPl  =  PI *2.000 

CONS T  9.020-14*1 IM*F)»»1 .5 

SNRHO  --DSINlRHO  ERAD) 

SNRHO  =  DSGR  T [ sN  RHO ) 

MSUM  ( 0 . ODO , 0 . ODO ) 

DO  45  MD*1,NM 
L  F  =  0 

23  INI  10 
LF  -=  LF  ♦  1 
DO  25  1  =  1  ,  Nf 

IF | MODE  1 MD , I  I  . E Q .  0)  GO  TO  25 
JJ  =  I-KKLMD)*1 
MF  =  MODE l MO, I) 

XXI  UJ  I  =  F Rl  0 1  1  ) 

YV I JJ  )  =  YCt  LF ,MD , J J ) 

Bl JJ 1  =  BCiLF.MD.JJ) 

Cl JJ )  =  CCl L I , MD , J J ) 

D ( J J 1  =  DC(LF.MD.JJ) 

25  CONTINUE 

IF  I F .GE . XX( 1 ) )  GO  TO  30 
GO  TO  45 

30  I F I F  .  LE . X X l M F  I ) GO  TO  33 
GO  TO  45 
33  CONTIN.E 

El LF  1 r  SPEVAL I F , XX . Y Y . B . C . D . MF ,  I N I T ) 

35  I  F  (  LF  .  LT  .4  )  L',0  T0  23 

CPXSIN=CDSIN( E(3)+IM»E(4) ) 

E  XC  - 1  I  1  I  ♦  ID-  L)S(  E(  2)  )  +  IM*DSIN(  E(2)  )  ) 

CAY  =  T  WOP  I  *  F  L  L 1 T  E 

MSUM-MSUM*E.\C*CDEXPl  -1M* CAY* RHO*  (CPXSIN-1  .  DO  )  ) 

45  CONTINUE 

CHN L= CONST *MSUM/ SNRHO 

RETURN 

END 
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FUNCTION  SPEVAL  (XVAL,  X,  Y  ,  0,  C.  0.  N.  I N I  T  ) 
IMPLICIT  REAL-Bt A-H.O'Z' 

DIMENSION  Til),  Y(1),  0(1).  C(1).  D(  1  ) 

SP  EVAL  EVALUA1ES  THE  INTERPOLATING  CUBIC  SPLINE 
FOR  THE  DATA  (Ml).Yll)).  1*1. ...N  AT  X  =  XVAL. 

IT  IS  ASSUMED  THAT  THE  CUBIC  POLYNOMIALS  GIVEN 
IN  Bill.  Cil).  DID  HAVE  BEEN  PREVIOUSLY 
COMPUTED  BY  THE  SUBROUTINE  SPLINE  OR  PSFL1N. 

1NIT  IS  AN  ESTIMATE  OF  ThE  INTERVAL  WHERE  XVAL 
LIES,  Xi INI T ). LE . AVAL . LE . X( INIT+1 ) ,  BUT  NEED 
NOT  BE  USED.  SET  1N1T=0  IF  THERE  IS  NO  ESTIMATE. 
ON  RETURN.  I N I  I  WILL  CONTAIN  THE  INTERVAL  NUMBER. 


FN 

=  N  -  1 

EPS 

■  *  1 . 0E-3  »  t  X  i  N  )  - 

Xi  1  )  ) 

/  fn 

IF 

l XV AL . LT . M 1  ) -E  PS ) 

GO  TO 

800 

IF 

l XV AL .GT . XiN l+EPS) 

GO  TO 

BOO 

IF 

i 1 N1 T . LE . 0 )  GO  TO 

200 

IF 

i INI T .GE. N)  GO  TO 

200 

I  F 

iXVAL.LT.XilNIT)) 

GO  TO 

150 

100  IF 

(XVAL. LT. XI INIT  +  1  )  ) 

GO  TO  300 

IF  i INIT+1 .GL.N)  GO  TO  300 


INI T  =  INIT  +  1 
GO  TO  TOO 

150  INIT  *  INIT  -  1 

IF  ( INIT . LE.O)  GO  TO  200 

IF  ( XV AL .GE. X( INI T ) )  GO  TO  300 

GO  TO  ISO 

200  INIT  =  1 
GO  TO  100 

300  H  *  XVAL  -  X i INIT) 

SPEVAL  =  Y( INIT )  ♦ 

$  (ID(INIT).H  ♦  Cl  INIT) ) *  H  ♦  B i  1 N I T ) ) *H 

RETJRN 

BOO  PRINT  BOO 

PRINT  901.  XVAL,X( 1 ) ,X(N) 

RETURN 

900  FORMAT  i'  ERROR  IN  SP  EVAL',/, 

$  '  XVAL  OUT  OF  INTERPOLATION  RANGE') 

901  FORMAT  15X,'  XVAL  * '  , 1  PE  1 2 . 5 ,  '  X(1)  »',1PE12.5.'  X(N)*',1PE12 

END 
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SUBROUTINE  NlOON  ( N , X , V , S I GNT , A , B ) 
IMPLICIT  REAl*61A-M.0-2) 

DIMENSION  X(  t) ,  Yl  1  ) ,  MI  15) 

LX  »  2  •  •  N 

F  L  X  LX 

FLXI 1 » ( B-A )  TLX 

P12  =6  .  283 1  BE' J07 

PI  24  *  I -PI  2  4.0  ) 

DO  1  I  =  1  ,  N 

I  Mill  =  2 • *  I N -  1) 

DO  4  L  >  1  ,N 
NULOCK  »  2**IL-1) 

LBLOCK  *  LX  \ 3 LOCK 
LBMALF  =  LB.OCK  2 
KO  =  0 

DO  4  I  BLOCK  1.  NBLOCK 
ISTART  »  LB.  l'CK«  l  1  BLOCK  -  1) 

F  ft  i  ftO 

V  =  ( S I GNT •  >’  l  2  •  FK  ) , FLX 
21 » DCOS I  V 1 
22*DSlNiV) 

IF  IL'ABSIV  -  P 1 24  )  -  1.0D-6)  11,  1 

II  22  -  -1.0 

12  CONI [Nut 

DO  2  I  •  1 .  LBMALF 
J  *  I  START  ♦  1 

K  *  J  *  LBMALF 
01  X i K  1  * 2 1  -  ' I K ) *22 
02  »  Y I K 1  * 2 1  ♦  X  I K ) • 22 
X(  K  1  3  XI  J  )  -  O' 

Y ( K 1  =  Y I J )  -  02 

X  i  U 1  =  X ( J  )  *  O' 

Y ( J 1  =  Y I U )  ♦  02 

2  CONTINUE 

DO  3  I  *  2,  N 

I  I  =  I 

LL  =  ANDl  M l I > , KO ) 

I F I L  L  )  4,4,3 

3  KO  ’  KO  -  Mil) 

4  KO  *  KO  ♦  Ml  I  I ) 

KO  =  0 

DO  50  K  *  1 .  LX 
K 1  =K0  +  1 

IF  i K 1 -K ) 55 , 55 , 65 
65  Ml  a  X  I  K 1  ) 

H2  =  Y I K 1 ) 

X  I  K  1)  *  X  I  K  ) 

Y  I  K  1  )  »  Y  l  K  ) 

XlKl^HI 

Y I K ) 3  H2 

55  CONTINUE 

DO  65  I  »  1  ,  N 

II  ^  I 
LL  ■ 


AND! M( I ) , KO ) 


IF  Ul  )  75 

,7b, 85 

85 

KO  = 

KO  - 

Mil) 

75 

KO  - 

kO  ♦ 

Mill) 

50 

COM 

INUE 

DO  10 

0  K  =  1 

.  1  X 

OV  =  1'I 

M  K- 

11  t  B- 

A) 

l  F  l  K 

.  NE  . 

1  1  GO 

TO  9 

G  -  F  l  x 

1  1 

GO  10 

1  0 

9 

G 1 1  -1 

.  mDSINiOM* 

F  l  X  1  1 

G»G 

F  l  X  I  1 

10 

COM 

INUE 

X  t  K  1  = 

Xt  K  )  ♦ 

G 

j; 

) ) • »2 ) / l OM»OM ) 

’ 
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